---
title: "Evaluating the effect of irradiation on the densities of two RNA viruses in                                                 Glossina morsitans morsitans"
author: "Mirieri et al.,"
date: "13-7-2022"
output:
  word_document: default
  html_document: default
---

```{r}


setwd("C:/Users/abdallaa/OneDrive - IAEA/My_passport_6/Caroline_Wageningen/Viruses_manuscript/15022023/Iflanegemanuscript202317thFeb")

library(ggplot2)
library(MASS)
library(rmarkdown)
library(knitr)
library(lme4)
library(MuMIn)
library(ggthemes) # Load
library(datasets)
library(plyr)
library(dplyr)
library(tidyverse)

```



```{r}
#------------------------------------------------------

##Normality for neg, b and C

#Normality for neg

qpcr <- read.csv("Iflaneg_may_2023.csv",sep=",", row.names=NULL)

#qpcr$Life_stage=as.factor(qpcr$Life_stage)
qpcr$Irradiation_dose=as.factor(qpcr$Irradiation_dose)
str(qpcr)

attach(qpcr)
#head(neg)
qpcr=na.omit(qpcr)
qpcr

# check for normality of the data
Box=boxcox(Normalized_expression~1,  lambda = seq(-35,15,0.1))
Cox = data.frame(Box$x, Box$y)
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),]
Cox2[1,]
lambda = Cox2[1, "Box.x"] 
lambda

Transformed = (Normalized_expression ^ lambda - 1)/lambda
Transformed

expression <- with(qpcr,Normalized_expression)
trans <- with(qpcr,Transformed)
library(rcompanion)

plotNormalHistogram(Normalized_expression)
qqnorm(Normalized_expression,ylab="Normalized_expression")
qqline(Normalized_expression, 
       col="red")

plotNormalHistogram(Transformed)
qqnorm(Transformed,ylab="Transneg_expression")
qqline(Transformed, 
       col="red")
```



```{r}
#---------------------------------------------------------------

#===============================================================================

ifla<-subset(qpcr, Target=="GmmIV")


#select the model

anova(with(ifla,lm(Transformed ~ Irradiation_dose)))
anova(with(ifla,lm(Transformed ~ Life_stage)))
anova(with(ifla,lm(Transformed ~ Irradiation_conditions)))

anova(with(ifla,lm(Transformed ~ Irradiation_dose +Life_stage)))
anova(with(ifla,lm(Transformed ~ Irradiation_dose+Irradiation_conditions)))
anova(with(ifla,lm(Transformed ~ Irradiation_dose *Life_stage)))
anova(with(ifla,lm(Transformed ~ Irradiation_dose *Irradiation_conditions)))

#---------------------------------------------------------------------------------------------------------------


mod1 <- lm(Transformed ~ Irradiation_dose + Irradiation_conditions, data = ifla)
mod2 <- lm(Transformed ~ Irradiation_dose * Irradiation_conditions, data = ifla)
mod3 <- lm(Transformed ~ Irradiation_dose + Life_stage + Irradiation_conditions, data = ifla)
mod4 <- lm(Transformed ~ Irradiation_dose * Life_stage * Irradiation_conditions, data = ifla)
mod5 <- lm(Transformed ~ Irradiation_dose + Life_stage, data = ifla)
mod6 <- lm(Transformed ~ Irradiation_dose*Life_stage, data = ifla)
mod7 <- lm(Transformed ~ Irradiation_dose, data = ifla)
mod8 <- lm(Transformed ~ Irradiation_conditions, data = ifla)
mod9 <- lm(Transformed ~  Life_stage, data = ifla)
mod10<- lm(Transformed ~ Irradiation_dose*Life_stage, data = ifla)

AICc(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9,mod10)
anova(mod8, mod9)

##no difference, choose either model but went for the lowest

summary(mod8)

#------------------------------------------------------------------------------




```




```{r}



#====================================================================================
###Irradiation dose

ifla$Irradiation_dose=as.factor(ifla$Irradiation_dose)


Fig1a<-ggplot(ifla,aes(x=Irradiation_dose,y=Transformed, fill=Irradiation_dose)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)
Fig1a

tiff("Fig 1a.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(Fig1a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_blank()) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation doses (Gy)"))) + ylab(expression (bold("Transformed relative  GmmIV density")))
dev.off()


###Irradiation status
#ifla$Irradiation_conditions=as.factor(ifla$Irradiation_conditions)

Fig1b<-ggplot(ifla,aes(x=Irradiation_conditions,y=Transformed, fill=Irradiation_conditions)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)
Fig1b

tiff("Fig 1b.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(Fig1b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_blank()) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation status "))) + ylab(expression (bold("Transformed  relative GmmIV density")))
dev.off()


###Days_post.irradiation
#ifla$Life_stage=as.factor(ifla$Life_stage)

Fig1c<-ggplot(ifla,aes(x=Life_stage,y=Transformed, fill=Life_stage)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)
Fig1c

tiff("Fig 1c.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(Fig1c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_blank()) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Days post irradiation"))) + ylab(expression (bold("Transformed  relative GmmIV density")))
dev.off()



###Irradiation dose according to air status 


Fig2a<-ggplot(ifla, aes(x=Irradiation_dose,y=Transformed, colour=Irradiation_conditions)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

Fig2a


tiff("Fig 2a.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(Fig2a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose"))) + ylab(expression(bolditalic("GmmIV_ ")~bold("Transformed normalized density")))
dev.off()



###Irradiation dose according to pupae stage/ days post irradiation
Fig2b<-ggplot(ifla, aes(x=Irradiation_dose,y=Transformed, colour=Life_stage)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

Fig2b


tiff("Fig 2b.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(Fig2b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose(Gy)"))) + ylab(expression(bolditalic("GmmIV_ ")~bold("Transformed normalized density")))
dev.off()


```


```{r}

# #Statistics iflaevirus  for subsets 


#-------------------------------------------------------------------

#iflaevirus-Normoxia
#---------------------------------------------
#iflaNorm<- read.csv("iflaevirus_alldata_normalized.csv")
iflaNorm<- subset(ifla, Irradiation_conditions=="Normoxia")

#-------------------------------------------------------------------
anova(with(iflaNorm,lm(Transformed~ Life_stage*Irradiation_dose)))
anova(with(iflaNorm,lm(Transformed~ Life_stage +Irradiation_dose)))
anova(with(iflaNorm,lm(Transformed~ Life_stage)))
anova(with(iflaNorm,lm(Transformed~ Irradiation_dose)))


mod1 <- lm(Transformed ~ Irradiation_dose, data = iflaNorm)
mod5 <- lm(Transformed~ Irradiation_dose+Life_stage, data = iflaNorm)
mod6 <- lm(Transformed ~ Irradiation_dose*Life_stage, data = iflaNorm)
mod11 <- lm(Transformed ~Life_stage, data = iflaNorm)


AICc(mod1, mod5, mod6, mod11)  
anova(mod11, mod1)

summary(mod11)



#-------------------------------------------------------------

#FIGURE_normoxia
#iflaNorm$Life_stage=as.factor(iflaNorm$Life_stage)

FigS2A<-ggplot(iflaNorm, aes(x=Irradiation_dose,y=Transformed, colour= Life_stage
)) +geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

FigS2A


tiff("Fig S2A.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(FigS2A+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose"))) + ylab(expression(bolditalic("GmmIV_ ")~bold("Normoxia transformed relative density")))
dev.off()




#------------------------------------------------------
#iflaevirus-Hypoxia    
#-----------------------------------------------
#iflaHypox <- read.csv("iflaevirus_alldata_normalized.csv")
iflaHypox <- subset(ifla , Irradiation_conditions=="Hypoxia")

#-------------------------------------------------------------------------
anova(with(iflaHypox,lm(Transformed~  Life_stage+Irradiation_dose)))
anova(with(iflaHypox,lm(Transformed~  Life_stage*Irradiation_dose)))
anova(with(iflaHypox,lm(Transformed~  Life_stage)))
anova(with(iflaHypox,lm(Transformed~  Irradiation_dose)))


mod1 <- lm(Transformed ~ Irradiation_dose, data = iflaHypox)
mod5 <- lm(Transformed~ Irradiation_dose+Life_stage, data = iflaHypox)
mod6 <- lm(Transformed ~ Irradiation_dose*Life_stage, data = iflaHypox)
mod11 <- lm(Transformed ~Life_stage, data = iflaHypox)

AICc(mod1, mod5, mod6, mod11)
anova(mod11, mod6)
summary(mod11)
summary(mod6)

#HSD analysis in case of significant differenes
iflahypoxdata1<-aov(with(data=iflaHypox,lm(Transformed ~ Irradiation_dose*Life_stage)))
iflahypoxdata2 <- TukeyHSD(iflahypoxdata1,ordered = FALSE, conf.level = 0.95)
iflahypoxdata2

#---------------------------------------------------------------

###Figure hypoxia
FigS2B<-ggplot(iflaHypox, aes(x=Irradiation_dose,y=Transformed, colour= Life_stage)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

FigS2B


tiff("Fig S2B.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(FigS2B+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose (Gy)"))) + ylab(expression(bolditalic("GmmIV_ ")~bold("Hypoxia Transformed relative density")))
dev.off()






```





```{r}
#
#work with GmmNegV data

neg<-subset(qpcr, Target=="GmmNegV")

#neg$Life_stage=as.factor(neg$Life_stage)
neg$Irradiation_dose=as.factor(neg$Irradiation_dose)
#select the model

anova(with(neg,lm(Transformed ~ Irradiation_dose)))
anova(with(neg,lm(Transformed ~ Life_stage)))
anova(with(neg,lm(Transformed ~ Irradiation_conditions)))

anova(with(neg,lm(Transformed ~ Irradiation_dose + Irradiation_conditions)))
anova(with(neg,lm(Transformed ~ Irradiation_dose+Life_stage)))
anova(with(neg,lm(Transformed ~ Irradiation_dose*Irradiation_conditions)))
anova(with(neg,lm(Transformed ~ Irradiation_dose*Life_stage)))
#---------------------------------------------------------------------------------------------------------------


mod1 <- lm(Transformed ~ Irradiation_dose + Irradiation_conditions, data = neg)
mod2 <- lm(Transformed ~ Irradiation_dose * Irradiation_conditions, data = neg)
mod3 <- lm(Transformed ~ Irradiation_dose + Life_stage + Irradiation_conditions, data = neg)
mod4 <- lm(Transformed ~ Irradiation_dose * Life_stage * Irradiation_conditions, data = neg)
mod5 <- lm(Transformed ~ Irradiation_dose + Life_stage, data = neg)
mod6 <- lm(Transformed ~ Irradiation_dose*Life_stage, data = neg)
mod7 <- lm(Transformed ~ Irradiation_dose, data = neg)
mod8 <- lm(Transformed ~ Irradiation_conditions, data = neg)
mod9 <- lm(Transformed ~  Life_stage, data = neg)
mod10<- lm(Transformed ~ Irradiation_dose+Life_stage, data = neg)

AICc(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9, mod10)
anova(mod8, mod9)

##no difference, choose either model but went for the lowest

summary(mod8)
summary(mod9)
#------------------------------------------------------------------------------

#====================================================================================

neg$Irradiation_dose=as.factor(neg$Irradiation_dose)


Fig3a<-ggplot(neg,aes(x=Irradiation_dose,y=Transformed, fill=Irradiation_dose)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)
Fig3a

tiff("Fig 3a.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(Fig3a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_blank()) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation doses (Gy)"))) + ylab(expression (bold("Transformed relative  GmmNegV density")))
dev.off()


###Irradiation_conditions
#neg$Irradiation_conditions=as.factor(neg$Irradiation_conditions)

Fig3b<-ggplot(neg,aes(x=Irradiation_conditions,y=Transformed, fill=Irradiation_conditions)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)
Fig3b

tiff("Fig 3b.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(Fig3b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_blank()) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Irradiation status"))) + ylab(expression (bold("Transformed  relative GmmNegV density")))
dev.off()


###Days_post.irradiation
###neg$Life_stage=as.factor(neg$Life_stage)

Fig3c<-ggplot(neg,aes(x=Life_stage,y=Transformed, fill=Life_stage)) +
  geom_boxplot() + geom_jitter(width=0.1,alpha=0.2)
Fig3c

tiff("Fig 3c.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(Fig3c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(axis.text.x = element_blank()) +theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + xlab(expression(bold("Days post irradiation"))) + ylab(expression (bold("Transformed  relative GmmNegV density")))
dev.off()


###Irradiation dose according to irradiation status 


Fig4a<-ggplot(neg, aes(x=Irradiation_dose,y=Transformed, colour=Irradiation_conditions)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

Fig4a


tiff("Fig 4a.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(Fig4a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose(Gy)"))) + ylab(expression(bolditalic("GmmNegV_ ")~bold("Transformed normalized density")))
dev.off()




###Irradiation dose according to days post irradiation
Fig4b<-ggplot(neg, aes(x=Irradiation_dose,y=Transformed, colour=Life_stage)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

Fig4b


tiff("Fig 4b.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(Fig4b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose(Gy)"))) + ylab(expression(bolditalic("GmmNegV_ ")~bold("Transformed normalized density")))
dev.off()



```




```{r}


#-------------------------------------------------------------------

#GmmNegV-Normoxia
#---------------------------------------------
#NegNorm<- read.csv("GmmNegV_alldata_normalized.csv")
NegNorm<- subset(neg, Irradiation_conditions=="Normoxia")

#-------------------------------------------------------------------
anova(with(NegNorm,lm(Transformed~ Life_stage*Irradiation_dose)))
anova(with(NegNorm,lm(Transformed~ Life_stage +Irradiation_dose)))
anova(with(NegNorm,lm(Transformed~  Life_stage)))
anova(with(NegNorm,lm(Transformed~  Irradiation_dose)))

mod1 <- lm(Transformed ~ Irradiation_dose, data = NegNorm)
mod5 <- lm(Transformed~ Irradiation_dose+Life_stage, data = NegNorm)
mod6 <- lm(Transformed ~ Irradiation_dose*Life_stage, data = NegNorm)
mod11 <- lm(Transformed ~Life_stage, data = NegNorm)


AICc(mod1, mod5, mod6, mod11) 
anova(mod1, mod5)

summary(mod1)

#HSD analysis in case of significant differenes
data1<-aov(with(data=NegNorm,lm(Transformed~  Irradiation_dose)))
data2 <- TukeyHSD(data1,ordered = FALSE, conf.level = 0.95)
data2

#-------------------------------------------------------------

#FIGURE_normoxia
###NegNorm$Life_stage=as.factor(NegNorm$Life_stage)

FigsS3A<-ggplot(NegNorm, aes(x=Irradiation_dose,y=Transformed, colour= Life_stage
)) +geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

FigsS3A


tiff("Fig S3A.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(FigsS3A+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose (Gy)"))) + ylab(expression(bolditalic("GmmNegV_ ")~bold("Normoxia transformed relative density")))
dev.off()

#-----------------------------------------------
#NegHypox <- read.csv("GmmNegV_alldata_normalized.csv")
NegHypox <- subset(neg , Irradiation_conditions=="Hypoxia")

#-------------------------------------------------------------------------
anova(with(NegHypox,lm(Transformed~  Life_stage+Irradiation_dose)))
anova(with(NegHypox,lm(Transformed~  Life_stage*Irradiation_dose)))
anova(with(NegHypox,lm(Transformed~  Life_stage)))
anova(with(NegHypox,lm(Transformed~  Irradiation_dose)))


mod1 <- lm(Transformed ~ Irradiation_dose, data = NegHypox)
mod5 <- lm(Transformed~ Irradiation_dose+Life_stage, data = NegHypox)
mod6 <- lm(Transformed ~ Irradiation_dose*Life_stage, data = NegHypox)
mod11 <- lm(Transformed ~Life_stage, data = NegHypox)

AICc(mod1, mod5, mod6, mod11)

anova(mod11,mod5)
summary(mod11)




#---------------------------------------------------------------

###Figure hypoxia
FigS3B<-ggplot(NegHypox, aes(x=Irradiation_dose,y=Transformed, colour= Life_stage)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

FigS3B


tiff("Fig S3B.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(FigS3B+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title= element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Irradiation dose (Gy)"))) + ylab(expression(bolditalic("GmmNegV_ ")~bold("Hypoxia Transformed relative density")))
dev.off()








```






